This document contains all analyses on the hanging chain behavior in weaver ants (Oecophylla smaragdina). Raw data manipulation not included.

setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

if (!require(librarian)) {
  install.packages("librarian")
  library(librarian)
}
## Loading required package: librarian
shelf(tidyverse, lme4, survival, glmmTMB, ggeffects, sjPlot, lmerTest, deSolve, DHARMa, survminer, data.table, viridis, scales, patchwork, profvis, COEF, ggh4x)
## 
##   The 'cran_repo' argument in shelf() was not set, so it will use
##   cran_repo = 'https://cran.r-project.org' by default.
## 
##   To avoid this message, set the 'cran_repo' argument to a CRAN
##   mirror URL (see https://cran.r-project.org/mirrors.html) or set
##   'quiet = TRUE'.
## Warning in checkDepPackageVersion(dep_pkg = "TMB"): Package version inconsistency detected.
## glmmTMB was built with TMB version 1.9.1
## Current TMB version is 1.9.4
## Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)
# Load predict method for tls plots
predictdf.TLS <- function(model, xseq, se, level) {
  pred <- as.numeric(model[[1]] %*% t(cbind(1, xseq)))
  if(se) {
    preds <- sapply(1:length(model[[2]]), function(x) as.numeric(c(model[[2]][x],model[[3]][x]) %*% t(cbind(1, xseq))))
    data.frame(
      x = xseq,
      y = pred,
      ymin = apply(preds, 1, function(x) quantile(x, probs = (1-level)/2)),
      ymax = apply(preds, 1, function(x) quantile(x, probs = 1-((1-level)/2)))
    )
  } else {
    return(data.frame(x = xseq, y = pred))
  }
}

. MAIN ANALYSES

. Position of joining over the chain

## Load dataset
join_dist_reg <-read.csv("joining_distance.csv")

## Calculate proportions of attaching ants at the various points of the chain (area under curve)
DensityFunction = approxfun(density(join_dist_reg$tip_distance))

AreaUnderCurve = function(lower, upper) {
    integrate(DensityFunction, lower=lower, upper=upper) }

## Find 95% of data
quantile(join_dist_reg$tip_distance, probs = c(0.025, 0.975))
##      2.5%     97.5% 
## -5.014442 14.977483
x <- join_dist_reg$tip_distance

MinToZero <- round(with(density(x, from = -7, to = 0), sum(y * diff(x)[1])), digits = 2)*100
ZeroToTen <- round(with(density(x, from = 0, to = 10), sum(y * diff(x)[1])), digits = 2)*100
TenToMax <- round(with(density(x, from = 10), sum(y * diff(x)[1])), digits = 2)*100

redline <- density(x, from = -5.014442, to = 14.977483)

a <- density(x, from = round(min(x)), to = 0)
b <- density(x, from = 0, to = 10)
c <- density(x, from = 10)

ChainWalked_plot_upd <- ggplot() +
  geom_line(data.frame(a), mapping = aes(x = a$x, y = a$y)) +
  geom_line(data.frame(b), mapping = aes(x = b$x, y = b$y)) +
  geom_line(data.frame(c)[c$x <= max(x),], mapping = aes(x = x, y = y)) +
  geom_line(data.frame(redline), mapping = aes(x = redline$x, y = redline$y), col = "red", lwd = 1) +
  geom_area(aes(x = stage(join_dist_reg$tip_distance,
                          after_stat = scales::oob_censor(x, c(-7, 0))),
                fill = "#0000FFA0"),
            stat = "density",
            alpha = .4) +
  theme_pubr(base_size = 20,
             legend = "none") + 
    geom_area(aes(x = stage(join_dist_reg$tip_distance,               
                          after_stat = scales::oob_censor(x, c(-0, 10))),
                fill = "#A0A0A0A0"),
            stat = "density",
            alpha = .4) +
    geom_area(aes(x = stage(join_dist_reg$tip_distance, 
                          after_stat = scales::oob_censor(x, c(10, max(join_dist_reg$tip_distance)))),
                fill = "#FF0000A0"),
            stat = "density",
            alpha = .4) +
  xlab("Distance from chain tip (mm)") +
  geom_segment(aes(x = mean(join_dist_reg$tip_distance), y = 0, xend = mean(join_dist_reg$tip_distance), yend = 0.11), 
               color = "blue", linetype = 2, size = 1) +
    labs(y = "density") +
  scale_fill_brewer(palette = "Dark2",                 ## color blind vision friendly
                    name = "Percentage of ants",
                    direction = -1,
                    guide = "legend",
                    labels = c(paste0(MinToZero, "%"), paste0(ZeroToTen, "%"), paste0(TenToMax, "%"))) + 
  guides(fill = guide_legend(override.aes = list(alpha = 0.2))) +
  theme_pubr(base_size = 20,
             legend = c(0.75, 0.84)) +
    scale_x_continuous(breaks = c(round(min(join_dist_reg$tip_distance)), 0, 10, 20),
                     labels = c(round(min(join_dist_reg$tip_distance)), 0, 10, 20)) +
  geom_segment(aes(x = -5.014442, y = 0, xend = 14.977483, yend = 0), lwd = 1, col = "red") +
  geom_segment(aes(x = -5.014442, y = 0, xend = -5.014442, yend = 0.0337), lwd = 1, col = "red") +
  geom_segment(aes(x = 14.977483, y = 0, xend = 14.977483, yend = 0.011), lwd = 1, col = "red") 


ChainWalked_plot_upd

. Probability of joining beyond the tip of the chain

## Load dataset
joining_distance <- read.csv("joining_distance.csv") %>%
  mutate(scaled.GapHeight = scale(GapHeight, scale = T, center = T), # center variables
         scaled.gaplength = scale(gaplength, scale = T, center = T),
         lengthen = ifelse(Scaled_y > tip_y, 1, 0)) # binarize response

# Fit logistic model
lengthen_mod <- glmmTMB(lengthen ~ scaled.GapHeight + scaled.gaplength + scaled.GapHeight : scaled.gaplength + (1|RepName), data = joining_distance, family = binomial)

# Model diagnostic
res <- simulateResiduals(lengthen_mod, n = 1000)
testQuantiles(res)

## 
##  Test for location of quantiles via qgam
## 
## data:  simulationOutput
## p-value = 0.9049
## alternative hypothesis: both
testDispersion(res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.005, p-value = 0.906
## alternative hypothesis: two.sided
testZeroInflation(res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0034, p-value = 1
## alternative hypothesis: two.sided
testUniformity(res)

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.065447, p-value = 0.7806
## alternative hypothesis: two-sided
testOutliers(res)

## 
##  DHARMa bootstrapped outlier test
## 
## data:  res
## outliers at both margin(s) = 0, observations = 96, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0 0
## sample estimates:
## outlier frequency (expected: 0 ) 
##                                0
# See model results
summary(lengthen_mod)
##  Family: binomial  ( logit )
## Formula:          
## lengthen ~ scaled.GapHeight + scaled.gaplength + scaled.GapHeight:scaled.gaplength +  
##     (1 | RepName)
## Data: joining_distance
## 
##      AIC      BIC   logLik deviance df.resid 
##    130.4    143.2    -60.2    120.4       91 
## 
## Random effects:
## 
## Conditional model:
##  Groups  Name        Variance  Std.Dev. 
##  RepName (Intercept) 6.488e-09 8.055e-05
## Number of obs: 96, groups:  RepName, 15
## 
## Conditional model:
##                                   Estimate Std. Error z value Pr(>|z|)  
## (Intercept)                        -0.2710     0.2560  -1.059   0.2898  
## scaled.GapHeight                   -0.1619     0.3007  -0.538   0.5902  
## scaled.gaplength                   -0.5433     0.3090  -1.758   0.0787 .
## scaled.GapHeight:scaled.gaplength  -0.2686     0.3202  -0.839   0.4016  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Get model predictions
pred_lengthen_gl <- ggeffect(lengthen_mod, terms = "scaled.gaplength") 

## Unscale variable before plotting
pred_lengthen_gl$unscaled <- pred_lengthen_gl$x * attr(joining_distance$scaled.gaplength, 'scaled:scale') + attr(joining_distance$scaled.gaplength, 'scaled:center')

# Create plot object
lengthen_plot <- ggplot(pred_lengthen_gl, aes(unscaled, predicted)) +
  geom_line(size = .7) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1, fill = "red") +
  theme_classic() +
  labs(x = "Distance from target platform (mm)", y = "probability to lengthen the chain") +
  ylim(0, 1) +
  theme_pubr(margin = T, base_size = 20, base_family = "Times New Roman") 

## Visualize plot
lengthen_plot

. Analysis of leaving positions of ants

leaving_positions <- read.csv("leaving_positions.csv")

leavedist <- leaving_positions$distancefromtip

al <- density(leavedist, to = 10)
bl <- density(leavedist, from = 10)

TipArea <- round(with(al, sum(y * diff(x)[1])), digits = 2)*100
OtherArea <- round(with(bl, sum(y * diff(x)[1])), digits = 2)*100
redline <- density(leavedist, from = 0, to = 10)

leave_position_plot <- ggplot() +
  geom_density(data = leaving_positions, aes(x = distancefromtip)) +
  geom_area(aes(x = stage(leaving_positions$distancefromtip,
                          after_stat = scales::oob_censor(x, c(10, 30))),
                fill = "#0000FFA0"),
            stat = "density",
            alpha = .4) +
  geom_area(aes(x = stage(leaving_positions$distancefromtip,               
                          after_stat = scales::oob_censor(x, c(0, 10))),
                fill = "#A0A0A0A0"),
            stat = "density",
            alpha = .4) +
  geom_segment(aes(x = mean(leavedist), y = 0, xend = mean(leavedist), yend = 0.165), 
               color = "blue", linetype = 2, size = 1) +
  labs(x = "Distance from chain tip (mm)", 
       y = "density") +
  scale_fill_brewer(palette = "Dark2",                 ## color blind vision friendly
                    name = "Percentage of ants",
                    direction = -1,
                    guide = "legend",
                    labels = c(paste0(OtherArea, "%"), paste0(TipArea, "%"))) + 
  guides(fill = guide_legend(override.aes = list(alpha = 0.2,
                                                 color = NA),
                             reverse = T)) +
  theme_pubr(base_size = 20,
             legend = c(0.75, 0.88)) +
    scale_x_continuous(breaks = c(round(min(leavedist)), 0, 10, 20),
                     labels = c(round(min(leavedist)), 0, 10, 20),
                     limits = c(0, 25)) ; leave_position_plot
## Warning: Removed 205 rows containing missing values (`position_stack()`).
## Warning: Removed 307 rows containing missing values (`position_stack()`).

## Time spent by ants that we followed
timespent <- read.csv("timespent.csv")

summary(timespent$duration)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   3.276   6.025  11.229  23.198  30.251 117.833
mean(timespent$duration)
## [1] 23.19815
sd(timespent$duration)
## [1] 26.34429
min(timespent$duration)
## [1] 3.276
max(timespent$duration)
## [1] 117.833

. Analysis of chain length change before and after leaving event (from elsewhere than the tip)

chain_growth <- read.csv("chain_growth.csv")
leave_10 <- read.csv("leave_10.csv")

chain_length_L <- data.frame()
for (i in 1:nrow(leave_10)) {
  
  df <- leave_10[i,]
  before <- mean(chain_growth[chain_growth$Rep == df$Rep & 
                      between(chain_growth$Slice, df$from, df$Slice),]$funct_y)*-1
  after <- mean(chain_growth[chain_growth$Rep == df$Rep & 
                      between(chain_growth$Slice, df$Slice, df$to),]$funct_y)*-1
  
  out <- data.frame(ID = df$ID,
                    before = before, 
                    after = after)
  
  chain_length_L <- rbind(chain_length_L, out)
  
  
}

chain_length_w <- chain_length_L %>%
  mutate(change = after - before) %>%
  pivot_longer(names_to = "when", cols = c(before, after)) %>%
  mutate(funct_y = value) 

ggplot(chain_length_w, aes(x = when, y = funct_y, group = ID, color = ID)) +
  geom_pointpath() +
  scale_x_discrete(limits=rev, labels = c("Before leaving event",
                                          "After leaving event"),
                   expand = expansion(mult = 0.2)) +
  ylim(15, 50) +
  theme_538(base_size = 15) +  
  theme(legend.position = "none",
        panel.border = element_rect(colour = "black", fill=NA, linewidth=2)) +
  labs(x = "",
       y = "Length of chain (mm)") +
  geom_text(data = chain_length_w[chain_length_w$when == "before" & 
                                    chain_length_w$ID != "A7_32_50" & 
                                    chain_length_w$ID != "A3_50_50",], 
            aes(y = funct_y + 0.3, x = 2.07, label = round(change, digits = 2))) +
  geom_text(data = chain_length_w[chain_length_w$when == "before" &
                                    chain_length_w$ID == "A7_32_50",], 
            aes(y = funct_y + 0.5, x = 2.07, label = round(change, digits = 2))) + 
  geom_text(data = chain_length_w[chain_length_w$when == "before" &
                                    chain_length_w$ID == "A3_50_50",], 
            aes(y = funct_y - 0.5, x = 2.07, label = round(change, digits = 2)))

mean(chain_length_w[chain_length_w$when == "before",]$change)
## [1] 0.213505
sd(chain_length_w[chain_length_w$when == "before",]$change)
## [1] 0.2706782

. Probability of joining and leaving the structure

. Joining probability

# Load dataset
join_growth <- read.csv("JL_chaingrowth.csv") %>%
  filter(Behavior != "Extension") %>%
  mutate(scaled.gaplength = scale(gaplength_j, center = T, scale = T),
         scaled.GapHeight = scale(GapHeight, center = T, scale = T))

# Fit logistic model
join_mod <- glmer(Join ~ scaled.gaplength + scaled.GapHeight + scaled.gaplength:scaled.GapHeight + traf_join_down + traf_join_up + (1|RepName), 
                  data = join_growth,
                  family = binomial,
                  glmerControl(optimizer = "Nelder_Mead"))

# Residual diagostics
res <- simulateResiduals(join_mod, n = 1000)
testUniformity(res)

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.058243, p-value = 0.2247
## alternative hypothesis: two-sided
testDispersion(res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.0031, p-value = 0.928
## alternative hypothesis: two.sided
testQuantiles(res)

## 
##  Test for location of quantiles via qgam
## 
## data:  simulationOutput
## p-value = 0.3492
## alternative hypothesis: both
testZeroInflation(res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0026, p-value = 1
## alternative hypothesis: two.sided
# Results
summary(join_mod)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## Join ~ scaled.gaplength + scaled.GapHeight + scaled.gaplength:scaled.GapHeight +  
##     traf_join_down + traf_join_up + (1 | RepName)
##    Data: join_growth
## Control: glmerControl(optimizer = "Nelder_Mead")
## 
##      AIC      BIC   logLik deviance df.resid 
##    451.4    477.8   -218.7    437.4      315 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3896 -1.1302  0.7864  0.8516  1.4849 
## 
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  RepName (Intercept) 0.02371  0.154   
## Number of obs: 322, groups:  RepName, 15
## 
## Fixed effects:
##                                   Estimate Std. Error z value Pr(>|z|)  
## (Intercept)                        0.39051    0.17654   2.212    0.027 *
## scaled.gaplength                   0.16773    0.15961   1.051    0.293  
## scaled.GapHeight                  -0.03996    0.17682  -0.226    0.821  
## traf_join_down                     0.15286    0.53416   0.286    0.775  
## traf_join_up                      -0.82711    0.78701  -1.051    0.293  
## scaled.gaplength:scaled.GapHeight -0.24377    0.18630  -1.308    0.191  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) scld.g scl.GH trf_jn_d trf_jn_p
## scld.gplngt  0.324                                
## scld.GpHght -0.170 -0.621                         
## traf_jn_dwn -0.422 -0.001 -0.199                  
## traf_join_p -0.241  0.008 -0.178  0.130           
## scld.gp:.GH -0.428 -0.597  0.643 -0.099   -0.126
# Plot 
## Predict values from model
pred_join_gl <- ggeffect(join_mod, terms = "scaled.gaplength")

## Unscale variables before plotting
pred_join_gl$unscaled <- pred_join_gl$x * attr(join_growth$scaled.gaplength, 'scaled:scale') + attr(join_growth$scaled.gaplength, 'scaled:center')

# Create plot for gaplength
join_plot1 <- ggplot(pred_join_gl, aes(unscaled, predicted)) +
  geom_line(size = .5) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1, fill = "red") +
  ylim(0, 1) +
  theme_classic() +
  labs(x = "Distance from platform (mm)", y = "Probability to join the chain") 

# Forest plot of the model
join_plot2 <- plot_model(join_mod,
                         title = "",
                         show.p = TRUE,
                         axis.lim=c(.08, 5), 
                         dot.size = 4, 
                         line.size = 1.6) +
  scale_x_discrete(labels=c(
    scaled.GapHeight = "Experimental condition",
    scaled.gaplength = "Distance from platform",
    'scaled.gaplength:scaled.GapHeight' = "Experimental condition*\nDistance from platform",
    traf_join_down = "Traffic down",
    traf_join_up = "Traffic up")) +
  scale_color_brewer(palette = "Dark2", direction = -1) +
  theme_pubr(base_size = 15) +
  scale_y_continuous(labels = c(0.1, 1, 2, 3), 
                     breaks = c(0.1, 1, 2, 3)) +
  geom_hline(aes(yintercept = 1), linetype = 2);join_plot2
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

## Visualize plot
join_plot1

join_plot2

. Leaving probability

## Load dataset
leave_df <- read.csv("JL_chaingrowth.csv") %>%
  filter(Behavior != "Extension", 
         Join == 1)  # Only ants that joined can leave

## Calculate number of ants joining while focal ant is in the structure
njoin_df <- list()
for (i in 1:nrow(leave_df)) {
  df <- leave_df[i,]
  Rep <- df$RepName
  start <- df$tstart_l
  stop  <- df$tstop_l
  df$njoin <- nrow(leave_df[between(leave_df$tstop_j, start, stop) & 
                              leave_df$Subject != df$Subject & 
                              leave_df$RepName == Rep,])
  njoin_df[[length(njoin_df) + 1]] <- df
}

## Re-create updated dataframe
leave_df <- do.call(rbind, njoin_df)

## Tidy dataframe
leave_df <- leave_df %>%
  select(c("GapHeight", "RepName", "Subject", "Leave", "time_spent", "traf_leave_down", "traf_leave_up", "gaplength_l", "njoin"))

## Transform njoin in rate (Rjoin) to normalize it over time
leave_df$Rjoin <- leave_df$njoin/leave_df$time_spent

## Create surv object
leave_surv <- Surv(leave_df$time_spent, leave_df$Leave, type = "right")

## Fit model
leave_mod <- coxph(Surv(leave_df$time_spent, leave_df$Leave, type = "right") ~ gaplength_l + traf_leave_down + traf_leave_up + Rjoin, data = leave_df)

## COx model diagnostics
## PH assumptions are respected
cox.zph(leave_mod)
##                 chisq df     p
## gaplength_l     0.236  1 0.627
## traf_leave_down 2.411  1 0.121
## traf_leave_up   1.510  1 0.219
## Rjoin           0.401  1 0.527
## GLOBAL          9.154  4 0.057
## Results
summary(leave_mod)
## Call:
## coxph(formula = Surv(leave_df$time_spent, leave_df$Leave, type = "right") ~ 
##     gaplength_l + traf_leave_down + traf_leave_up + Rjoin, data = leave_df)
## 
##   n= 180, number of events= 42 
## 
##                       coef  exp(coef)   se(coef)      z Pr(>|z|)    
## gaplength_l      6.529e-02  1.067e+00  1.901e-02  3.436 0.000591 ***
## traf_leave_down -1.047e+01  2.831e-05  2.479e+00 -4.225 2.39e-05 ***
## traf_leave_up    5.542e+00  2.551e+02  1.469e+00  3.772 0.000162 ***
## Rjoin           -8.957e-01  4.083e-01  3.230e+00 -0.277 0.781514    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                 exp(coef) exp(-coef) lower .95 upper .95
## gaplength_l     1.067e+00  9.368e-01 1.028e+00 1.108e+00
## traf_leave_down 2.831e-05  3.532e+04 2.199e-07 3.645e-03
## traf_leave_up   2.551e+02  3.919e-03 1.433e+01 4.543e+03
## Rjoin           4.083e-01  2.449e+00 7.277e-04 2.291e+02
## 
## Concordance= 0.819  (se = 0.036 )
## Likelihood ratio test= 62.86  on 4 df,   p=7e-13
## Wald test            = 70.54  on 4 df,   p=2e-14
## Score (logrank) test = 82.97  on 4 df,   p=<2e-16
## Visualize effects
## Effect of gaplength on leaving probability
newdata_gl <- data.frame(gaplength_l = c(quantile(leave_df$gaplength_l)[[1]], 
                                         quantile(leave_df$gaplength_l)[[5]]),
                         traf_leave_up = mean(leave_df$traf_leave_up), 
                         traf_leave_down   = mean(leave_df$traf_leave_down),
                         Rjoin = mean(leave_df$Rjoin))

# plot
leaveplot_gl <- ggsurvplot(survfit(leave_mod, newdata = newdata_gl), 
                           data = leave_df,
                           conf.int = F,
                           title = "Distance from target platform",
                           legend.title = "",
                           legend.labs = c("5.29mm",
                                           "39.29mm"),
                           ylim = c(0, 1),
                           xlim = c(0, 75),
                           palette = c("#d95f02", "#7570b3"),
                           size = 1.5)[[1]] +
  theme_pubr(base_size = 20,
             legend = c(0.5, 0.1)) +
  theme(legend.key.size = unit(.7, "cm"),
        plot.title = element_text(hjust = 0.5, face = "bold"),
        legend.text = element_text(size = 20),
        legend.title = element_blank(),
        legend.direction = "horizontal"); leaveplot_gl
## Warning: `gather_()` was deprecated in tidyr 1.2.0.
## ℹ Please use `gather()` instead.
## ℹ The deprecated feature was likely used in the survminer package.
##   Please report the issue at <https://github.com/kassambara/survminer/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Effect of traffic up on leaving probability
newdata_tfup <- data.frame(gaplength_l = mean(leave_df$gaplength_l),
                           traf_leave_up = c(quantile(leave_df$traf_leave_up)[[1]],
                                             quantile(leave_df$traf_leave_up)[[5]]), 
                           traf_leave_down   = mean(leave_df$traf_leave_down),
                           Rjoin = mean(leave_df$njoin))

# plot
leaveplot_tfup <- ggsurvplot(survfit(leave_mod, newdata = newdata_tfup), 
                             data = leave_df, 
                             conf.int = F,
                             title = "Traffic leaving chain's tip",
                             legend.title = "Leaving traffic rate:",
                             legend.labs = c("0 ants/s",
                                             "0.57 ants/s"),
                             ylim = c(0, 1),
                             xlim = c(0, 75),
                             palette = c("#d95f02", "#7570b3"),
                             size = 1.5,
                             legend = c(0.5, 0.05))[[1]]  +
  theme_pubr(base_size = 20,
             legend = c(0.5, 0.1)) +
  theme(legend.key.size = unit(.7, "cm"),
        plot.title = element_text(hjust = 0.5, face = "bold"),
        legend.text = element_text(size = 20),
        legend.title = element_blank(),
        legend.direction = "horizontal"); leaveplot_tfup

# Effect of traffic down on leaving probability
newdata_tfdw <- data.frame(gaplength_l = mean(leave_df$gaplength_l),
                           traf_leave_up   = mean(leave_df$traf_leave_up),
                           traf_leave_down = c(quantile(leave_df$traf_leave_down)[[1]],
                                               quantile(leave_df$traf_leave_down)[[5]]),
                           Rjoin = mean(leave_df$njoin))

# plot
leaveplot_tfdw <- ggsurvplot(survfit(leave_mod, newdata = newdata_tfdw), data = leave_df, 
                          conf.int = F,
                          title = "Traffic arriving at chain's tip",
                          legend.title = "Arriving traffic rate:",
                          legend.labs = c("0 ants/s",
                                          "0.59 ants/s"),
                          xlim = c(0, 75),
                          palette = c("#d95f02", "#7570b3"),
                          size = 1.5,
                          legend = c(0.5, 0.05))[[1]] +
  theme_pubr(base_size = 20,
             legend = c(0.5, 0.1)) +
  theme(legend.key.size = unit(.7, "cm"),
        plot.title = element_text(hjust = 0.5, face = "bold"),
        legend.text = element_text(size = 20),
        legend.title = element_blank(),
        legend.direction = "horizontal"); leaveplot_tfdw

. Reaching behaviour

# Load dataset
Extension_df <- read.csv("JL_chaingrowth.csv") %>%
  filter(Behavior == "Extension") %>%
  filter(tstop_j - tstart_j != 0) %>% 
  mutate(scaled.GapHeight = scale(GapHeight, center = T, scale = T),
         scaled.gaplength = scale(gaplength_j, center = T, scale = T))

# Fit logistic model
extend_mod <- glmer(Extend ~ scaled.GapHeight + scaled.gaplength + scaled.gaplength : scaled.GapHeight + traf_join_up + traf_join_down + (1|RepName), 
                    data = Extension_df, 
                    family = binomial, 
                    glmerControl(optimizer = "Nelder_Mead"))

## Residual diagnostics
res <- simulateResiduals(extend_mod, n = 1000)
testUniformity(res)

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.045241, p-value = 0.515
## alternative hypothesis: two-sided
testDispersion(res)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.9996, p-value = 0.986
## alternative hypothesis: two.sided
testQuantiles(res)

## 
##  Test for location of quantiles via qgam
## 
## data:  simulationOutput
## p-value = 0.2627
## alternative hypothesis: both
testZeroInflation(res)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0013, p-value = 1
## alternative hypothesis: two.sided
## See results
summary(extend_mod)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## Extend ~ scaled.GapHeight + scaled.gaplength + scaled.gaplength:scaled.GapHeight +  
##     traf_join_up + traf_join_down + (1 | RepName)
##    Data: Extension_df
## Control: glmerControl(optimizer = "Nelder_Mead")
## 
##      AIC      BIC   logLik deviance df.resid 
##    368.2    394.8   -177.1    354.2      320 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.9672 -0.6154 -0.4699  1.1062  2.9815 
## 
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  RepName (Intercept) 0.1095   0.3309  
## Number of obs: 327, groups:  RepName, 15
## 
## Fixed effects:
##                                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                       -1.07573    0.21982  -4.894  9.9e-07 ***
## scaled.GapHeight                   0.04711    0.21689   0.217  0.82803    
## scaled.gaplength                  -0.47992    0.18276  -2.626  0.00864 ** 
## traf_join_up                       0.20555    0.63210   0.325  0.74505    
## traf_join_down                    -0.36416    0.60118  -0.606  0.54468    
## scaled.GapHeight:scaled.gaplength -0.11688    0.21508  -0.543  0.58682    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) scl.GH scld.g trf_jn_p trf_jn_d
## scld.GpHght -0.104                                
## scld.gplngt  0.367 -0.509                         
## traf_join_p -0.263 -0.051 -0.063                  
## traf_jn_dwn -0.411 -0.103 -0.079  0.020           
## scld.GpHg:. -0.362  0.658 -0.507 -0.009   -0.031
# Prediction for plotting
pred_extend <- ggeffect(extend_mod, terms = c("scaled.gaplength"))

# Unscale variable before plotting
pred_extend$unscaled <- pred_extend$x * attr(Extension_df$scaled.gaplength, "scaled:scale") + attr(Extension_df$scaled.gaplength, "scaled:center")

extend_plot1 <- ggplot(pred_extend, aes(unscaled, predicted)) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1, fill = "red") +
  theme_pubr(base_size = 20) +
  labs(x = "Distance from target platform (mm)", y = "Probability of displaying reaching behavior") 

extend_plot1

. Latency to join the chain is not modulated by state of the structure

# Load dataset 
joinlatency_df <- read.csv("JL_chaingrowth.csv") %>%
  filter(Behavior != "Extension") %>%
  filter(Join == 1) 

# Create survival object
join_latency <- Surv(joinlatency_df$latency_tojoin, joinlatency_df$Join, type = "right")

# Fit survival model
join_latency_mod <- coxph(join_latency ~ gaplength_latency + GapHeight:gaplength_latency + traf_join_down + traf_join_up, data = joinlatency_df)

# Diagnostics
cox.zph(join_latency_mod)
##                              chisq df    p
## gaplength_latency           0.5356  1 0.46
## traf_join_down              0.0453  1 0.83
## traf_join_up                0.0308  1 0.86
## gaplength_latency:GapHeight 0.6362  1 0.43
## GLOBAL                      0.8086  4 0.94
# Results
summary(join_latency_mod)
## Call:
## coxph(formula = join_latency ~ gaplength_latency + GapHeight:gaplength_latency + 
##     traf_join_down + traf_join_up, data = joinlatency_df)
## 
##   n= 180, number of events= 180 
## 
##                                   coef  exp(coef)   se(coef)      z Pr(>|z|)
## gaplength_latency            0.0196222  1.0198160  0.0321975  0.609    0.542
## traf_join_down               0.0367758  1.0374604  0.4635245  0.079    0.937
## traf_join_up                 0.2993326  1.3489582  0.6486954  0.461    0.644
## gaplength_latency:GapHeight -0.0001222  0.9998778  0.0005739 -0.213    0.831
## 
##                             exp(coef) exp(-coef) lower .95 upper .95
## gaplength_latency              1.0198     0.9806    0.9574     1.086
## traf_join_down                 1.0375     0.9639    0.4182     2.574
## traf_join_up                   1.3490     0.7413    0.3783     4.810
## gaplength_latency:GapHeight    0.9999     1.0001    0.9988     1.001
## 
## Concordance= 0.529  (se = 0.026 )
## Likelihood ratio test= 2.51  on 4 df,   p=0.6
## Wald test            = 2.55  on 4 df,   p=0.6
## Score (logrank) test = 2.55  on 4 df,   p=0.6

. Analytical model

. Parameter estimation from experiments

. Probability of joining

The probability of joining is constant over time. We calculate the baseline probability of joining from the statistical model at the average conditions of the chain.

## Probability of joining the chain at average values
average_data <- data.frame(
  scaled.gaplength = mean(join_growth$scaled.gaplength),
  scaled.GapHeight = mean(join_growth$scaled.GapHeight),
  traf_join_down = mean(join_growth$traf_join_down),
  traf_join_up = mean(join_growth$traf_join_up))

## This is the baseline probability of joining the chain.
joining_probability <- predict(join_mod, newdata = average_data, re.form = ~0, type = "response")

joining_probability
##        1 
## 0.587615

. Probability of leaving the chain

The probability of leaving depends largely on the distance from the visual target.

# Calculate decay rate for each experiment
pred <- list()
for (GL in seq(0, max(leave_df$gaplength_l), 1)) {
  
  # New dummy dataframe with variable number of ants
  newdata_prediction <- data.frame(gaplength_l = GL,
                                   traf_leave_up = mean(leave_df$traf_leave_up), 
                                   traf_leave_down = mean(leave_df$traf_leave_down),
                                   Rjoin = mean(leave_df$Rjoin))
  
  
  # Build dataframe with survival time and survival probability
  survival_df <- data.frame(survfit(leave_mod, newdata = newdata_prediction, censor = T, start.time = 0)$time,
                            survfit(leave_mod, newdata = newdata_prediction, censor = T, start.time = 0)$surv)
  
  # rename columns
  colnames(survival_df) <- c("time", "surv") 
  
  survival_df <- survival_df %>%
    filter(time < 80)
  
  # Take log of survival probability
  survival_df$survlog <- log(survival_df$surv) 
  
  # Fit linear model to find constant decay rate
  survival_logmod <- lm(survlog ~ time, survival_df)

  # Extract slope
  slope <- survival_logmod$coefficients[["time"]]

  # Write in list
  pred[[length(pred) + 1]] <- data.frame(gaplength = GL, R = slope)

}

# Transform in dataframe
predictions <- do.call(rbind, pred)

# Take log of the survival estimates (linarization)
predictions$logR <- log(-predictions$R)

# Fit linear regression
R_logmod <- lm(logR ~ gaplength, predictions)

# Extract coefficients from the model
Pl0 <- R_logmod$coefficients[1]
slope <- R_logmod$coefficients[2]

# Fit exponential curve to model predictions
w_plot <- ggplot(predictions, aes(gaplength, -R)) +
  geom_point() +
  geom_line(aes(x = gaplength, y = exp(Pl0)*exp(slope*gaplength)), color = "blue", size = 1) +
  labs(x = "Distance from target platform (mm)",
       y = "W (decay rate)") +
  theme_pubr() 
  
w_plot

. Alternative model decay rate

In the alternative model, where the probability of leaving is independent of S(t), we calculate w using the estimated survival probability of ants at the average conditions of the chain.

# Average chain conditions
newdata_prediction <- data.frame(gaplength_l = mean(leave_df$gaplength_l),
                                 traf_leave_up = mean(leave_df$traf_leave_up),
                                 traf_leave_down = mean(leave_df$traf_leave_down),
                                 Rjoin = mean(leave_df$Rjoin))

# Build dataframe with survival time and survival probability
survival_df <- data.frame(survfit(leave_mod, newdata = newdata_prediction, censor = T, start.time = 0)$time,
                          survfit(leave_mod, newdata = newdata_prediction, censor = T, start.time = 0)$surv)
# rename columns
colnames(survival_df) <- c("time", "surv")

survival_df <- survival_df %>%
  filter(time < 80)

# Take log of survival probability
survival_df$survlog <- log(survival_df$surv)

# Fit linear model to find constant decay rate
survival_logmod <- lm(survlog ~ time, survival_df)

# Extract slope
slope.avg <- survival_logmod$coefficients[["time"]]

wplot.complete <- w_plot +
  geom_line(aes(y = slope.avg*-1), color = "dark green", size = 1, linetype = 1) +
  theme_pubr(base_size = 20); wplot.complete

. Relationship between number of ants and length of the chain

cumsum_df <- read.csv("cumsum_df.csv")

# Fit linear regression
mod_nlength <- lmer(chainlength ~ Nants + (1|RepName), cumsum_df)

## Model diagnostics
res <- simulateResiduals(mod_nlength, n = 1000)
testResiduals(res)

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.086222, p-value = 0.07049
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.0487, p-value = 0.72
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  simulationOutput
## outliers at both margin(s) = 1, observations = 225, p-value = 0.3624
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.0001125173 0.0245127474
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.004444444
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.086222, p-value = 0.07049
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.0487, p-value = 0.72
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  simulationOutput
## outliers at both margin(s) = 1, observations = 225, p-value = 0.3624
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.0001125173 0.0245127474
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.004444444
# Visualize fit lm
Nants_length_plot <- ggplot(cumsum_df, aes(Nants, chainlength)) +
  geom_point() +
  geom_line(aes(y = fixef(mod_nlength)[1] + fixef(mod_nlength)[2]*Nants), color = "blue", size = 1) +
  theme_pubr(base_size = 20) +
  xlab("Number of ants in the chain") +
  ylab("Length of the chain (mm)") 

Nants_length_plot

. Fit analytical model

Here we fit the analytical models.

# Load data on chain growth
growth_data <- read.csv("growth_data.csv") %>%
  separate("RepName", "GapHeight", sep = "_", remove = F, extra = "drop") %>%
  arrange(RepName, to)

# Loop 
chain.size = list()
for ( Rep in unique(growth_data$RepName)) {
  
  # Extract Experimental Condition from data
  G <- as.numeric(growth_data[growth_data$RepName == Rep,]$GapHeight[1])
  
  # Find maximum duration of each experiment
  duration = tail(ceiling(growth_data[growth_data$RepName == Rep,]$start), 1)
  
  # Extract traffic
   ## Traffic is used as rate over 10 seconds
   ## We input the exact traffic for each replicate, which will change every 10 seconds 
  Tr.Rep = growth_data[growth_data$RepName == Rep,] %>%
    select(start, to, traf_bottom_down) %>%
    mutate(traffic.rate = traf_bottom_down/10)
  
  # Probability of joining is fixed 
  Pj = joining_probability
  
  # Estimated coefficients for decay rate calculation
  R = R_logmod$coefficients[2]
  Pl.complete = exp(R_logmod$coefficients[1])

  # Set chain size at 0 (no ants in the chain)
  Size.t = 0 
  Size.t.avg = 0
  
  # Set time step of 0.5 seconds
  timestep = 0.5
  for (t in seq(10, duration, timestep)) {

    # Extract traffic rate 
    Tr = Tr.Rep[t >= Tr.Rep$start & t < Tr.Rep$to,]$traffic.rate
    
    # Calculate length of the chain for current value of S
    L = fixef(mod_nlength)[1] + fixef(mod_nlength)[2]*Size.t

    # Calculate distance from platform by removing estimated chain length from full gap 
    gaplength = G - L
    
    # Calculate decay rate for current gaplength
    w = Pl.complete*exp(R*gaplength)
    
    # Average model has fixed decay rate 
    w.avg = -slope.avg                                 
    
    # Find current size of chain
    # Our model
    Size = Size.t + (Pj * Tr - w * Size.t)*(timestep)

    # # Average model 
    Size.avg = Size.t.avg + (Pj * Tr - w.avg * Size.t.avg)*(timestep)

    # Store results for next time step
    Size.t = Size
    Size.t.avg = Size.avg
    
    # Export results
    chain.size[[length(chain.size) + 1]] <- data.frame(RepName = Rep,
                                                       GapHeight = G,
                                                       gaplength = (G-L),
                                                       chainlength_model = L,
                                                       time = t, 
                                                       S = Size,
                                                       S.avg = Size.avg)
    
  }
  
}

# List to dataframe
Size.avgmodel <- rbindlist(chain.size)

. Compare models

Compare.avgmodel <- merge(cumsum_df, Size.avgmodel, all.x = T, all.y = F, by = c("RepName", "time", "GapHeight")) %>%
  arrange(RepName, time) %>%
  # select(GapHeight, RepName, time, Nants, S, S.avg) %>%
  setnames("S", "model") %>%
  setnames("Nants", "experimental") %>%
  setnames("S.avg", "model.avg") %>%
  pivot_longer(cols = c(model, experimental, model.avg)) %>%
  setnames("name", "data") %>%
  setnames("value", "Nants")%>%
  mutate(labels = paste0("Target distance: ", GapHeight, "mm"),
         labels1 = paste0(GapHeight, "mm")) 

k = 1
dff <- list()
for (i in unique(Compare.avgmodel$RepName)){
  df <- Compare.avgmodel[Compare.avgmodel$RepName == i,]
  df$Rep <- paste0(df$GapHeight, "mm_Rep", k)
  dff[[length(dff) + 1]] <- df
  k <- k + 1
  if (k > 5) {k = 1}
}

Compare.avgmodel <- rbindlist(dff)

# Plot comparison
plotReps <- ggplot(Compare.avgmodel, aes(time, Nants, color = data)) +
  geom_point() +
  geom_line() +
  facet_wrap2("Rep", axes = T, remove_labels = F, ncol = 5) +
  labs(x = "Time (s)", 
       y = "Number of ants in the chain") +
  theme(legend.position = "right", 
        legend.title = element_blank()) +
  theme_pubr(base_size = 20) +
  scale_fill_brewer(palette = "Dark2", direction =-1,  aesthetics = c("colour", "fill"),
                    labels = c("Experimental", "Distance-dependent model", "Distance-independent model")) 

# Comparison between experimental and modelling data AVERAGE
avg.plot <- ggplot(Compare.avgmodel, aes(time, Nants, fill = data)) +
  stat_summary(fun.data = mean_se,geom = "point", aes(color = data)) +
  stat_summary(fun=mean,geom="line", aes(color = data)) +
  stat_summary(geom="ribbon", alpha = .4) +
  scale_fill_brewer(palette = "Dark2", direction = -1,  
                    aesthetics = c("colour", "fill"),
                    labels = c(" Experimental ", " Distance-dependent model ", " Distance-independent model "),
                    name = " ") +
  facet_wrap("labels1") +
  labs(x = "Time (s)", 
       y = "Number of ants in the chain") +
  theme(legend.title = element_blank()) +
  theme_classic() +
  theme_pubr(base_size = 20, legend = "top") +
  theme(strip.text.x = element_text(size = 20),
        legend.key.size = unit(.9, "cm"),
        legend.text = element_text(size = 20))  + 
  ylim(0,30) ; avg.plot
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`

# Create dataframe for estimatin glinear correlations
Compare.linear <- merge(cumsum_df, Size.avgmodel, all.x = T, all.y = F, by = c("RepName", "time", "GapHeight")) %>%
  arrange(RepName, time) %>%
  select(GapHeight, RepName, time, Nants, S, S.avg) %>%
  setnames("S", "model") %>%
  setnames("Nants", "experimental") %>%
  setnames("S.avg", "model.avg") %>%
  pivot_longer(cols = c(model, model.avg)) %>%
  setnames("name", "data") %>%
  setnames("value", "Nants") %>%
  mutate(labels = paste0("Target distance: ", GapHeight, "mm"),
         labels1 = paste0(GapHeight, "mm"))

linearplot <- ggplot(Compare.linear, aes(x = experimental, y = Nants, color = data, fill = data)) +
  geom_point() +
  geom_abline(aes(intercept = 0, slope = 1)) +
  facet_wrap("labels1") +
  theme_classic() +
  labs(title = "",
       x = "Experimental (n° of ants)",
       y     = "Model (n° of ants)") +
  geom_smooth(method = COEF::tls) +
  theme_pubr(base_size = 20, legend = "top") +
  scale_color_manual(values = c("#d95f02","#1b9e77"),
                     labels = c(" Distance-dependent model ", " Distance-independent model "),
                     name = " ") +
  scale_fill_manual(values = c("#d95f02","#1b9e77"),
                    labels = c(" Distance-dependent model ", " Distance-independent model "),
                    name = " ") + 
  theme(strip.text.x = element_text(size = 20), 
        legend.key.size = unit(.55, "cm"),
        legend.text = element_text(size = 20)) 
  
avg.plot
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`

linearplot 
## `geom_smooth()` using formula = 'y ~ x'

plotReps 

. Explore chains using modelling

# Assign names to variables
Pj = joining_probability
R = R_logmod$coefficients[2]
Pl.complete = exp(R_logmod$coefficients[1])
Tr.variance.avg = var(growth_data$traf_bottom_down/10)
max.time = max(growth_data$to)
slope.nlength <- fixef(mod_nlength)[2]
intercept.nlength <- fixef(mod_nlength)[1]

GapHeight = seq(10, 120, length.out = 41)
Traffic.mean = seq(0, 2, length.out = 41)
simulation.dataframe <- expand.grid(GapHeight, Traffic.mean)
colnames(simulation.dataframe) <- c("GapHeight", "Traffic.mean")
GapHeight.vec <- simulation.dataframe$GapHeight 
Traffic.mean.vec <- simulation.dataframe$Traffic.mean

# Function that outputs final chain length and time for each combination of traffic and GapHeight
simulateChain <- function(i, GH.vec, Tr.mean.vec, simulations) {
  
  GH = GH.vec[i]
  Tr.mean = Tr.mean.vec[i]
  
  # Calculate parameters for gamma distribution
  s = Tr.variance.avg/Tr.mean
  a = Tr.mean/s

  Tr.vec <- rgamma(
    n = simulations*501,
    shape = a,
    scale = s
  )
  # Counter for traffic
  counter.Tr = 1
  
  # Loop for simulations
  total.df <- vector("list", 1e5)
  for (sim in 1:simulations) {
    
    # Set chain size at 0 at the beginning of each simulation
    Size.t = 0
    
    # Start simulation
    rep.list <- vector("list", 1e5)
    # use counter to store values in list
    counter = 1
    time.step = 1
    
    for (t in seq(0, 500, time.step)) { ## 254s is the longest of my reps
      
      # Draw random traffic at each time step
      Tr <- Tr.vec[counter.Tr]
      
      # Calculate length of the chain
      Length.chain = intercept.nlength + slope.nlength*Size.t
      
      # Calculate distance from the platform
      gaplength = GH - Length.chain

      # Calculate decay rate 
      w = Pl.complete*exp(R*gaplength)
      w = ifelse(w > 1, 1, w)
      
      # Model
      Size = Size.t + (Pj * Tr - w * Size.t)*(time.step)
      
      # Save size at each time step
      Size.t <- Size
      if (Size.t < 0) {Size.t = 0}
      
      # increase counter
      counter = counter + 1
      counter.Tr = counter.Tr + 1

      # Store results of simulation
      rep.list[[counter]] <- list(time = t,
                                  S = Size.t,
                                  L = Length.chain,
                                  Tr = Tr.mean,
                                  GapHeight = GH)
      
      # Stop loop when chain length covers full gap
      if (Length.chain >= GH) {
        break
      }
      
    }
    
    # Output list 
    total.df[[sim]] <- rbindlist(rep.list)
    # Calculate slope for each simulation
    slope = coefficients(lm(L ~ time, total.df[[sim]]))[2]
    # Store slope results
    total.df[[sim]][[6]] <- slope
    names(total.df[[sim]])[[6]] <- "slope"
    # Store simulation number
    total.df[[sim]][[7]] <- sim
    names(total.df[[sim]])[[7]] <- "simulation"
    # Store iteration number
    total.df[[sim]][[8]] <- i
    names(total.df[[sim]])[[8]] <- "iteration"
    

  }
  
  # Put results in datatable
  tot.df <- rbindlist(total.df)
  
  # Output list
  tot.df
  
}

### Run model ### 
simulation.df <- rbindlist(lapply(1:length(GapHeight.vec), simulateChain, GH.vec = GapHeight.vec, Tr.mean.vec = Traffic.mean.vec, simulations = 100))

. Model results

## Function to calculate median growth rate
approxGrowth <- function(y, t) {
  median(diff(y[order(t)]))
}

dat <- simulation.df

dat <- dat %>%
  group_by(GapHeight, Tr, time) %>%
  summarise(L.avg = mean(L))
## `summarise()` has grouped output by 'GapHeight', 'Tr'. You can override using
## the `.groups` argument.
dat <- data.table(dat)

summ <- dat[, .(gr = approxGrowth(L.avg, time)), by = .(GapHeight, Tr)]

# Plot tile plot
SimulationTile_plot <- ggplot(summ[summ$Tr <= 1], aes(x = GapHeight, y = Tr)) +
  geom_tile(aes(fill = gr)) +
  stat_contour(aes(z = gr), color = "white", breaks = seq(0.01, 2, 0.25), linetype = 2) +
  stat_contour(aes(z = gr), color = "green", breaks = 0.01, linetype = 2) +
  coord_equal(ratio = max(summ$GapHeight) / max(summ$Tr)) +
  scale_x_continuous(expand=c(0,0)) +
  scale_y_continuous(expand=c(0,0)) +
  scale_fill_viridis_c(option = "inferno",
                       name = "growth 
rate (mm/s)") +
  labs(x = "Gap height (mm)", 
       y = "Traffic (ants/sec)",
       title = " ") +
  theme_minimal(base_size = 20) +
  theme(legend.title = element_text(),
        legend.justification = c(0, 1),
        legend.margin = margin(100, 0 , 0, 0),
        legend.spacing.y = unit(10, "pt")); SimulationTile_plot

ggpar(SimulationTile_plot,
      font.main = c(20, "bold"),
      font.x = c(20, "bold"),
      font.y = c(20, "bold"),
      font.caption = c(20, "bold"),
      
      font.tickslab = c(20, "bold"))

Myexps <- join_growth %>%
  group_by(RepName) %>%
  summarise(GapHeight = GapHeight[1],
    mean.Tr = mean(traf_join_down))

SimulationTile_plot +
 geom_point(data = Myexps, aes(GapHeight, mean.Tr), color = "white")

. Predictions from alternative model

# Assign names to variables
Pj = joining_probability
R = R_logmod$coefficients[2]
Pl.complete = exp(R_logmod$coefficients[1])
Tr.variance.avg = var(growth_data$traf_bottom_down/10)
max.time = max(growth_data$to)
slope.nlength <- fixef(mod_nlength)[2]
intercept.nlength <- fixef(mod_nlength)[1]
slope.avg <- slope.avg

# Create dataframe with values for simulation
GapHeight = seq(10, 120, length.out = 41)
Traffic.mean = seq(0, 2, length.out = 41)
simulation.dataframe <- expand.grid(GapHeight, Traffic.mean)
colnames(simulation.dataframe) <- c("GapHeight", "Traffic.mean")
GapHeight.vec <- simulation.dataframe$GapHeight 
Traffic.mean.vec <- simulation.dataframe$Traffic.mean

# Function that outputs final chain length and time for each combination of traffic and GapHeight
simulateChain <- function(i, GH.vec, Tr.mean.vec, simulations) {
  
  GH = GH.vec[i]
  Tr.mean = Tr.mean.vec[i]
  
  # Calculate parameters for gamma distribution
  s = Tr.variance.avg/Tr.mean
  a = Tr.mean/s

  Tr.vec <- rgamma(
    n = simulations*501,
    shape = a,
    scale = s
  )
  # Counter for traffic
  counter.Tr = 1
  
  # Loop for simulations
  total.df <- vector("list", 1e5)
  for (sim in 1:simulations) {
    
    # Set chain size at 0 at the beginning of each simulation
    Size.t = 0
    
    # Start simulation
    rep.list <- vector("list", 1e5)
    # use counter to store values in list
    counter = 1
    time.step = 1
    
    for (t in seq(0, 500, time.step)) { ## 254s is the longest of my reps
      
      # Draw random traffic at each time step
      Tr <- Tr.vec[counter.Tr]
      
      # Calculate length of the chain
      Length.chain = intercept.nlength + slope.nlength*Size.t
      
      # Calculate distance from the platform
      gaplength = GH - Length.chain
      
      # Calculate decay rate 
      w = -slope.avg ## change sign to obtain positive decay rate
      
      # Model
      Size = Size.t + (Pj * Tr - w * Size.t)*(time.step)
      
      # Save size at each time step
      # Size.t = ifelse(Size < 0, 0, Size)
      Size.t <- Size
      if (Size.t < 0) {Size.t = 0}
      
      # increase counter
      counter = counter + 1
      counter.Tr = counter.Tr + 1

      # Store results of simulation
      rep.list[[counter]] <- list(time = t,
                                  S = Size.t,
                                  L = Length.chain,
                                  Tr = Tr.mean,
                                  GapHeight = GH)
      
      # Stop loop when chain length covers full gap
      if (Length.chain >= GH) {
        break
      }
      
    }
    
    # Output list 
    total.df[[sim]] <- rbindlist(rep.list)
    # Calculate slope for each simulation
    slope = coefficients(lm(L ~ time, total.df[[sim]]))[2]
    # Store slope results
    total.df[[sim]][[6]] <- slope
    names(total.df[[sim]])[[6]] <- "slope"
    # Store simulation number
    total.df[[sim]][[7]] <- sim
    names(total.df[[sim]])[[7]] <- "simulation"
    # Store iteration number
    total.df[[sim]][[8]] <- i
    names(total.df[[sim]])[[8]] <- "iteration"
    

  }
  
  # Put results in datatable
  tot.df <- rbindlist(total.df)
  
  # Output list
  tot.df
  
}

### Run model ### 
simulation.avg.df <- rbindlist(lapply(1:length(GapHeight.vec), simulateChain, GH.vec = GapHeight.vec, Tr.mean.vec = Traffic.mean.vec, simulations = 100))
## Function to calculate median growth rate
approxGrowth <- function(y, t) {
  median(diff(y[order(t)]))
}

## Load simulated data
dat1<- simulation.avg.df

dat1 <- dat1 %>%
  group_by(GapHeight, Tr, time) %>%
  summarise(L.avg = mean(L))
## `summarise()` has grouped output by 'GapHeight', 'Tr'. You can override using
## the `.groups` argument.
dat1 <- data.table(dat1)

summ1 <- dat1[, .(gr = approxGrowth(L.avg, time)), by = .(GapHeight, Tr)]

# Plot tile plot
SimulationTile_avg_plot <- ggplot(subset(summ1, Tr <= 1)) +
  aes(x = GapHeight, y = Tr) +
  geom_tile(aes(fill = gr)) +
  stat_contour(aes(z = gr), color = "white", breaks = seq(0.01, 2, 0.25), linetype = 2) +
  stat_contour(aes(z = gr), color = "green", breaks = 0.01, linetype = 2) +
  coord_equal(ratio = max(summ1$GapHeight) / max(summ1$Tr)) +
  scale_x_continuous(expand=c(0,0)) +
  scale_y_continuous(expand=c(0,0)) +
  scale_fill_viridis_c(option = "inferno",
                      direction = 1,
                       name = "growth 
rate (mm/s)") +
  labs(x = "Gap height (mm)", 
       y = "Traffic (ants/sec)",
       title = " ") +
  theme_minimal(base_size = 16) +
  theme(legend.title = element_text(),
        legend.justification = c(0, 1),
        legend.margin = margin(100, 0 , 0, 0),
        legend.spacing.y = unit(10, "pt")) 

SimulationTile_avg_plot +
  geom_point(data = Myexps, aes(GapHeight, mean.Tr), color = "white")

. Time needed to complete a chain over a 110mm gap at the average traffic flow measured during our experiments

GapHeight.vec = 110
Traffic.mean.vec = 0.12

### Run model ### 
simulation.myexps <- rbindlist(lapply(1:length(GapHeight.vec), simulateChain, GH.vec = GapHeight.vec, Tr.mean.vec = Traffic.mean.vec, simulations = 100))

## Function to calculate median growth rate
approxGrowth <- function(y, t) {
  median(diff(y[order(t)]))
}

## Load simulated data
dat.myexps <-  simulation.myexps

dat.myexps <- dat.myexps %>%
  group_by(GapHeight, Tr, time) %>%
  summarise(L.avg = mean(L))
## `summarise()` has grouped output by 'GapHeight', 'Tr'. You can override using
## the `.groups` argument.
dat.myexps <- data.table(dat.myexps)

summ.myexps <- dat.myexps[, .(gr = approxGrowth(L.avg, time)), by = .(GapHeight, Tr)]

## Time needed to complete chain at estimated growth rate of alternative model
(110/summ.myexps$gr)/60
## [1] 31.44898

Additional simulations with distance from visual target constant

# Assign names to variables
Pj = joining_probability
R = R_logmod$coefficients[2]
Pl.complete = exp(R_logmod$coefficients[1])
Tr.variance.avg = var(growth_data$traf_bottom_down/10)
max.time = max(growth_data$to)
slope.nlength <- fixef(mod_nlength)[2]
intercept.nlength <- fixef(mod_nlength)[1]

# Create dataframe with values for simulation
GapHeight = seq(10, 120, length.out = 41)
Traffic.mean = seq(0, 2, length.out = 41)
simulation.dataframe <- expand.grid(GapHeight, Traffic.mean)
colnames(simulation.dataframe) <- c("GapHeight", "Traffic.mean")
GapHeight.vec <- simulation.dataframe$GapHeight 
Traffic.mean.vec <- simulation.dataframe$Traffic.mean
 
# Function that outputs final chain length and time for each combination of traffic and GapHeight
simulateChain <- function(i, GH.vec, Tr.mean.vec, simulations) {
  
  GH = GH.vec[i]
  Tr.mean = Tr.mean.vec[i]
  
  # Calculate parameters for gamma distribution
  s = Tr.variance.avg/Tr.mean
  a = Tr.mean/s

  Tr.vec <- rgamma(
    n = simulations*501,
    shape = a,
    scale = s
  )
  # Counter for traffic
  counter.Tr = 1
  
  # Loop for simulations
  total.df <- vector("list", 1e5)
  for (sim in 1:simulations) {
    
    # Set chain size at 0 at the beginning of each simulation
    Size.t = 0
    
    # Start simulation
    rep.list <- vector("list", 1e5)
    # use counter to store values in list
    counter = 1
    time.step = 1
    
    for (t in seq(0, 500, time.step)) { ## 254s is the longest of my reps
      
      # Draw random traffic at each time step
      Tr <- Tr.vec[counter.Tr]
      
      # Calculate length of the chain
      Length.chain = intercept.nlength + slope.nlength*Size.t
      
      # Calculate distance from the platform
      gaplength = 20
      
      # Calculate decay rate 
      w = Pl.complete*exp(R*gaplength)
      w = ifelse(w > 1, 1, w)
      
      # Model
      Size = Size.t + (Pj * Tr - w * Size.t)*(time.step)
      
      # Save size at each time step
      # Size.t = ifelse(Size < 0, 0, Size)
      Size.t <- Size
      if (Size.t < 0) {Size.t = 0}
      
      # increase counter
      counter = counter + 1
      counter.Tr = counter.Tr + 1

      # Store results of simulation
      rep.list[[counter]] <- list(time = t,
                                  S = Size.t,
                                  L = Length.chain,
                                  Tr = Tr.mean,
                                  GapHeight = GH)
      
      # Stop loop when chain length covers full gap
      if (Length.chain >= GH) {
        break
      }
      
    }
    
    # Output list 
    total.df[[sim]] <- rbindlist(rep.list)
    # Calculate slope for each simulation
    slope = coefficients(lm(L ~ time, total.df[[sim]]))[2]
    # Store slope results
    total.df[[sim]][[6]] <- slope
    names(total.df[[sim]])[[6]] <- "slope"
    # Store simulation number
    total.df[[sim]][[7]] <- sim
    names(total.df[[sim]])[[7]] <- "simulation"
    # Store iteration number
    total.df[[sim]][[8]] <- i
    names(total.df[[sim]])[[8]] <- "iteration"
    

  }
  
  # Put results in datatable
  tot.df <- rbindlist(total.df)
  
  # Output list
  tot.df
  
}

### Run model ### 
simulation.df_MP <- rbindlist(lapply(1:length(GapHeight.vec), simulateChain, GH.vec = GapHeight.vec, Tr.mean.vec = Traffic.mean.vec, simulations = 100))

dat_MovPl <- simulation.df_MP

dat_MovPl <- dat_MovPl %>%
  group_by(GapHeight, Tr, time) %>%
  summarise(L.avg = mean(L))
## `summarise()` has grouped output by 'GapHeight', 'Tr'. You can override using
## the `.groups` argument.
dat_MovPl <- data.table(dat_MovPl)

approxGrowth <- function(y, t) {
  median(diff(y[order(t)]))
}

summ_MovPl <- dat_MovPl[, .(gr = approxGrowth(L.avg, time)), by = .(GapHeight, Tr)]

# Plot tile plot
SimulationTile_MovPlatform <- ggplot(subset(summ_MovPl, Tr <= 1), aes(x = GapHeight, y = Tr)) +
  geom_tile(aes(fill = gr)) +
  stat_contour(aes(z = gr), color = "white", breaks = seq(0.01, 2, 0.25), linetype = 2) +
  stat_contour(aes(z = gr), color = "green", breaks = 0.01, linetype = 2) +
  coord_equal(ratio = max(summ_MovPl$GapHeight) / max(summ_MovPl$Tr)) +
  scale_x_continuous(expand=c(0,0)) +
  scale_y_continuous(expand=c(0,0)) +
  scale_fill_viridis_c(option = "inferno",
                       name = "growth 
rate (mm/s)") +
  labs(x = "Gap height (mm)", 
       y = "Traffic (ants/sec)",
       title = " ") +
  theme_minimal(base_size = 20) +
  theme(legend.title = element_text(),
        legend.justification = c(0, 1),
        legend.margin = margin(100, 0 , 0, 0),
        legend.spacing.y = unit(10, "pt")) 


ggpar(SimulationTile_MovPlatform,
      font.main = c(20, "bold"),
      font.x = c(20, "bold"),
      font.y = c(20, "bold"),
      font.caption = c(20, "bold"),
      
      font.tickslab = c(20, "bold"))

SimulationTile_MovPlatform

. Growth rate of chains when distance from platform is constant

GrowthMP <- read.csv("GrowthMP.csv") 
traffic_mean <- read.csv("traffic_mean.csv")

      ##### MODEL FITTING

Pj = joining_probability # Probability of joining the chain
duration <- max(GrowthMP$time_adj)
Tr.mean <- traffic_mean$traffic.mean
G.vec <- seq(2, 20, 1)
timestep = 1
Tr.variance.avg = traffic_mean$tr.var
s = Tr.variance.avg/Tr.mean  
a = Tr.mean/s

# Estimated coefficients for decay rate calculation
R = R_logmod$coefficients[2]
Pl.complete = exp(R_logmod$coefficients[1])

chain.size.MP = list()
for (sim in 1:1000) {
  
  Tr.vec <- rgamma( ## Traffic vector extracted from gamma distribution
    n = duration + 1,
    shape = a,
    scale = s
  )
  
  # Set chain size at 0 (no ants in the chain)
  Size.t = 0
  for (t in seq(1, duration, timestep)) {

    # Extract traffic rate
    Tr = Tr.vec[t]
    G = sample(G.vec, 1)

    # Calculate length of the chain for current value of S
    L = fixef(mod_nlength)[1] + fixef(mod_nlength)[2]*Size.t

    # Calculate decay rate for current gaplength
    w = Pl.complete*exp(R*G)
    
    # Find current size of chain
    # Our model
    Size = Size.t + (Pj * Tr - w * Size.t)*(timestep)
    
    # Store results for next time step
    Size.t = Size
    
    # Export results
    chain.size.MP[[length(chain.size.MP) + 1]] <- data.frame(time = t,
                                                       simulation = sim,
                                                       chainlength = L,
                                                       S = Size)
    
  }
  
}

# List to dataframe
model.MP <- rbindlist(chain.size.MP) 

model.MP <- model.MP %>%
  group_by(time) %>%
  mutate(
    low = quantile(chainlength, probs = 0.05, na.rm = T), 
    high = quantile(chainlength, probs = 0.95,  na.rm = T),
    type = "Model") 

GrowthMP <- GrowthMP %>%
  mutate(
    type = "Experimental"
  )

CompareMP_plot <- ggplot(model.MP, aes(x = time, y = chainlength, fill = type)) +
  stat_summary(data = GrowthMP, aes(time_adj, chainlength, col = type), fun.data = mean_se, inherit.aes = F) +
  stat_summary(fun.data = mean_se, geom = "line", size = 2, show.legend = T) +
  geom_ribbon(aes(ymin = low, ymax = high), alpha = 0.3, show.legend = T) +
  scale_fill_manual(name = "", values = c("#d95f02")) +
  scale_color_manual(name = " ", values = c("#000000")) +
  guides(color=guide_legend(override.aes=list(fill=NA, linetype=0)),
         fill=guide_legend(override.aes=list())) +
  xlab("Time (s)") +
  ylab("Length of the chain (mm)") +
  theme_pubr(base_size = 15) +
  labs(title = '(C) Constant-distance experiment') +
  theme(plot.title = element_text(face = "bold"),
        axis.title = element_text(size = 17),
        legend.key.size = unit(.9, "cm"),
        legend.text = element_text(size = 20),
        legend.position = c(0.1, 0.9),
        legend.spacing.y = unit(0, "cm")); CompareMP_plot
## Warning: Removed 8 rows containing missing values (`geom_segment()`).

. Plot arrangement and export

# Figure 2
leaveplot_tfdw_nl <- leaveplot_tfdw +
  theme(axis.title.y = element_blank(),
        axis.title.x = element_text(size = 20),
        title = element_text(size = 15)) 

leaveplot_gl_nl <- leaveplot_gl +
  theme(axis.title.y = element_blank(),
        axis.title.x = element_text(size = 20),
        title = element_text(size = 15))

leaveplot_tfup_nl <- leaveplot_tfup +
  theme(axis.title.x = element_text(size = 20),
        title = element_text(size = 15),
        axis.title.y = element_text(size = 15))

join_plot2_nl <- join_plot2 +
  theme(axis.title.x = element_text(size = 20), 
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 12)
        )

ChainWalked_plot_nl <- ChainWalked_plot_upd +
  ggtitle("Joining positions along chain") +
  theme(axis.text.x = element_text(size = 15),
        axis.title.x = element_text(size = 20),
        title = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        legend.title = element_text(size = 15),
        axis.title.y = element_text(size = 20),
        plot.title = element_text(face = "bold",
                                  hjust = 0.5)) 

## New plot for FIG 2

leave_position_plot_nl <- leave_position_plot +
  ggtitle("Leaving positions along chain") +
  theme(axis.text.x = element_text(size = 15),
        axis.title.x = element_text(size = 20),
        title = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        legend.title = element_text(size = 15),
        axis.title.y = element_text(size = 20),
        plot.title = element_text(face = "bold",
                                  hjust = 0.5)) 
  

Fig2 <- ggarrange(ggarrange(ChainWalked_plot_nl, leave_position_plot_nl,join_plot2_nl,ncol = 3, labels = c("(A)", "(B)", "(C)"), 
                            widths = c(0.35, 0.35, 0.3), hjust = c(-0.5, -0.5, -0.5)),
                  ggarrange(leaveplot_tfup_nl, leaveplot_tfdw_nl, leaveplot_gl_nl, ncol = 3, labels = ("(D)"), hjust = -0.5),
                  nrow = 2, heights = c(1, 0.85))

tiff("Fig2_upd2.tiff", units = "in", width = 15, height = 9, res = 900)
Fig2
dev.off()

tiff("Fig2_upd2.tiff", units = "in", width = 15, height = 9, res = 900)
Fig2
dev.off()

# Figure 3

avg.plot_nl <- avg.plot +
  scale_x_continuous(breaks = c(0, 50, 100, 150, 200, 250), labels = c(0, "",100,"", 200, "")) +
    theme(axis.title = element_text(size = 20), 
        legend.title = element_text(size = 0),
        legend.text = element_text(size = 20),
        legend.key.size = unit(.9, 'cm')) 

linearplot <- linearplot +
  scale_x_continuous(breaks = c(0, 5, 10, 15, 20, 25), labels = c(0, "",10,"", 20, "")) +
  theme(axis.title = element_text(size = 20), 
        legend.title = element_text(size = 0),
        legend.text = element_text(size = 20),
        legend.key.size = unit(.9, 'cm')) 


Fig3 <- ggarrange(avg.plot_nl, linearplot, common.legend = T, align = "v",
                  labels = c("(A)", "(B)"),font.label = list(size = 20),  nrow = 2) ; Fig3

tiff("Fig3.tiff", units = "in", width = 14, height = 20, res = 1200)
Fig3
dev.off()


# Figure 4
summ$Model <- "(A) Distance-dependent model"
summ1$Model <- "(B) Distance-independent model"
summ_MovPl$Model <- "(C) Constant-distance simulations"

CompareSim <- rbind(summ, summ1)
Myexps$Model <- "(A) Distance-dependent model"
Myexps1 <- Myexps %>%
  mutate(Model = "(B) Distance-independent model")
Myexps <- rbind(Myexps, Myexps1)

Fig4 <- ggplot(subset(CompareSim, Tr <= 1), aes(x = GapHeight, y = Tr)) +
  facet_wrap2("Model", ncol = 1, axes = T, remove_labels = F) +
  geom_tile(aes(fill = gr)) +
  stat_contour(aes(z = gr), color = "white", breaks = seq(0.25, 2, 0.25), linetype = 2) +
  stat_contour(aes(z = gr), color = "green", breaks = 0.01, linetype = 2) +
  scale_x_continuous(expand=c(0,0)) +
  scale_y_continuous(expand=c(0,0)) +
  scale_fill_viridis_c(option = "inferno",
                       name = "growth 
rate 
(mm/s)") +
  labs(x = "Gap height (mm)", 
       y = "Traffic (ants/sec)",
       title = " ") +
  theme_minimal(base_size = 12) +
  theme(legend.title = element_text(size = 12),
        legend.justification = c(0, 1),
        legend.margin = margin(0, 0 , 0, 0),
        strip.text = element_text(face = "bold", hjust = 0, size = 12))+
  geom_point(data = Myexps, aes(GapHeight, mean.Tr), color = "white")


CompareMP_plot <- ggplot(model.MP, aes(x = time, y = chainlength, fill = type)) +
  stat_summary(data = GrowthMP, aes(time_adj, chainlength, col = type), fun.data = mean_se, inherit.aes = F, size = 0.3) +
  stat_summary(fun.data = mean_se, geom = "line", size = 1, show.legend = T) +
  geom_ribbon(aes(ymin = low, ymax = high), alpha = 0.3, show.legend = T) +
  scale_fill_manual(name = "", values = c("#d95f02")) +
  scale_color_manual(name = " ", values = c("#000000")) +
  guides(color=guide_legend(override.aes=list(fill=NA, linetype=0)),
         fill=guide_legend(override.aes=list())) +
  xlab("Time (s)") +
  ylab("Length of the chain (mm)") +
  theme_pubr(base_size = 12) +
  labs(title = '  (C) Constant-distance experiment') +
  theme(plot.title = element_text(face = "bold", size = 12),
        axis.title = element_text(size = 12),
        legend.key.size = unit(.55, "cm"),
        legend.text = element_text(size = 10),
        legend.position = c(0.12, 0.9),
        legend.spacing.y = unit(0, "cm")); CompareMP_plot

Fig4_complete <- cowplot::plot_grid(Fig4, CompareMP_plot, 
                       ncol = 1, rel_heights = c(2, 1),
                       align = 'v', axis = 'lrtb') 


tiff("Fig4_MP.tiff", units = "cm", width = 18, height = 30, res = 600)
Fig4_complete
dev.off()

### SUPPLEMENTARY INFORMATION 

#Lengthening a chain (Fig S1)
tiff("FigS1.tiff", units = "in", width = 12, height = 10, res = 300)
lengthen_plot
dev.off()

# Reaching behavior

tiff("FigS2.tiff", units = "in", width = 12, height = 10, res = 300)
extend_plot1
dev.off()

# Log Linear model
FigS3 <- ggarrange(b, b1, b2, b3, labels = c("(A)", "(B)", "(C)", "(D)"),
          hjust = -0.1, vjust = 1.2)

tiff("FigS3.tiff", units = "in", width = 12, height = 12, res = 300)
FigS3
dev.off()


# W by distance from visual target

tiff("FigS4.tiff", units = "in", width = 12, height = 12, res = 300)
wplot.complete
dev.off()

# Length of chain vs. Number of ants 
tiff("FigS5.tiff", units = "in", width = 12, height = 12, res = 300)
Nants_length_plot
dev.off()

# Replicate plot v. modelling
tiff("FigS6.tiff", units = "in", width = 16, height = 12, res = 300)
plotReps
dev.off()

# Complete simulations plots
FigS8 <- ggplot(CompareSim, aes(x = GapHeight, y = Tr)) +
  facet_wrap2("Model", ncol = 2, axes = T, remove_labels = F) +
  geom_point(data = Myexps, aes(GapHeight, mean.Tr), color = "white") +
  geom_tile(aes(fill = gr)) +
  stat_contour(aes(z = gr), color = "white", breaks = seq(0.25, 2, 0.25), linetype = 2) +
  stat_contour(aes(z = gr), color = "green", breaks = 0.01, linetype = 2) +
  coord_equal(ratio = max(CompareSim$GapHeight) / max(CompareSim$Tr)) +
  scale_x_continuous(expand=c(0,0)) +
  scale_y_continuous(expand=c(0,0)) +
  scale_fill_viridis_c(option = "inferno",
                       name = "growth 
rate (mm/s)") +
  labs(x = "Gap height (mm)", 
       y = "Traffic (ants/sec)",
       title = " ") +
  theme_minimal(base_size = 15) +
  theme(legend.title = element_text(),
        legend.justification = c(0, 1),
        legend.margin = margin(0, 0 , 0, 0),
        strip.text = element_text(face = "bold", hjust = 0)) +
  geom_point(data = Myexps, aes(GapHeight, mean.Tr), color = "white")

  
tiff("FigS8.tiff", units = "in", width = 10, height = 6, res = 300)
FigS8
dev.off()